require(foreign)
require(dplyr)
require(stargazer)
require(reshape2)
require(MASS)
require(sandwich)

##All code run using R version 3.5.1##

###REPLICATION CODE FOR TRACKING TOLERANCE CHANGES####

datacfr <- read.csv("CFR_replication_data.csv")

##creating new variable for change in tolerance
datacfr$change <-  as.numeric(as.character(datacfr$EPA_2015_tolerance)) - 
  as.numeric(as.character(datacfr$EPA_1996_tolerance))

##variable for pesticide age
datacfr$year_new<-2015-datacfr$year_registered

#creating new factor variable for outcome of ordered logit
## 1 for less strict, 2 same, 3 stricter, 4 revoked (exempt means a tolerance no longer enforced)
datacfr$factor <- datacfr$change
datacfr$factor[datacfr$factor>0 |datacfr$EPA_2015_tolerance =="exempt"] <- 1
datacfr$factor[datacfr$factor==0] <- 2
datacfr$factor[datacfr$factor<0] <- 3
datacfr$factor[is.na(datacfr$factor)] <-4

factor <- datacfr$factor
datacfr$factor <- factor(factor)


###SUMMARY STATS FOR TABLE 1####

#Across all tolerances#
summary(datacfr$factor)/length(datacfr$factor)
length(datacfr$factor)
#For fruit/veggie
fv <- subset(datacfr,datacfr$Type=="fruit"|datacfr$Type=="vegetable")
summary(fv$factor)/length(fv$factor)
length(fv$factor)
#For grain
grain <-subset(datacfr,datacfr$Type=="grain")
summary(grain$factor)/length(grain$factor)
length(grain$factor)
#for animal product
animal <-subset(datacfr,datacfr$Type=="animal_product")
summary(animal$factor)/length(animal$factor)
length(animal$factor)
#pre 1961
pre_61<-subset(datacfr,datacfr$year_registered<1961)
summary(pre_61$factor)/length(pre_61$factor)
length(pre_61$factor)
#61-80
bet61_80<-subset(datacfr,datacfr$year_registered>=1961&datacfr$year_registered<1981)
summary(bet61_80$factor)/length(bet61_80$factor)
length(bet61_80$factor)
#post 1980
post_80<-subset(datacfr,datacfr$year_registered>1980)
summary(post_80$factor)/length(post_80$factor)
length(post_80$factor)
#most tox
most_tox<-subset(datacfr,datacfr$toxicity==5|datacfr$toxicity==4)
summary(most_tox$factor)/length(most_tox$factor)
length(most_tox$factor)
#least tox
least_tox<-subset(datacfr,datacfr$toxicity==1|datacfr$toxicity==2)
summary(least_tox$factor)/length(least_tox$factor)
length(least_tox$factor)


##MODELS FOR TABLE 2##
datacfr1<-subset(datacfr,!is.na(factor) & !is.na(year_new))
model_1 <- polr(factor~ year_new, data=datacfr1,Hess=TRUE)

##dropping NAs so model runs#
datacfrfull<-subset(datacfr1,!is.na(toxicity_worse) 
                    & !is.na(cancer_epa) 
                    & !is.na(toxicity) 
                    & !is.na(aquatic_chronic))
model_2<-polr(factor~ year_new + 
                toxicity+ 
                toxicity_worse+
                cancer_epa +
                aquatic_chronic +
                aquatic_acute, data=datacfrfull, Hess=TRUE)

model_full <-polr(factor~ year_new +
                    toxicity+ 
                    toxicity_worse+
                    cancer_epa +
                    aquatic_chronic +
                    aquatic_acute+
                    fruit_veg +
                    one_of_3_primary_acreage_crops_US + 
                    one_of_3_primary_acreage_crops_US*year_new
                  , Hess=TRUE,data= datacfrfull)


###models with robust standard errors###
vcov_1 <- vcovCL(model_1, as.character(datacfr1$pesticide))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_2, as.character(datacfrfull$pesticide))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_full, as.character(datacfrfull$pesticide))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_1,model_2,model_full, se=list(new.ses1, new.ses2, new.sesfull))



##ROBUSTNESS CHECK 1 (TABLE 3) - ELIMINATING DANGEROUS PESTICIDES####

data_safe_cfr<-subset(datacfr, datacfr$pesticide!="amitraz" & 
                        datacfr$pesticide!="azinphos-methyl" & 
                        datacfr$pesticide!="cyhexatin" & 
                        datacfr$pesticide!="disulfoton" & 
                        datacfr$pesticide!="endosulfan" & 
                        datacfr$pesticide!= "fenamiphos" & 
                        datacfr$pesticide!="isofenphos" & 
                        datacfr$pesticide!="methamidophos" & 
                        datacfr$pesticide!= "methidathion" & 
                        datacfr$pesticide!="mevinphos" & 
                        datacfr$pesticide!="phorate" & 
                        datacfr$pesticide!="terbufos")

######models for paper#####
data_safe_cfr<-subset(data_safe_cfr,!is.na(factor) & !is.na(year_new))
model_1 <- polr(factor~ year_new, data=data_safe_cfr)

data_safe_cfrfull<-subset(data_safe_cfr,!is.na(toxicity_worse) 
                          & !is.na(cancer_epa) 
                          & !is.na(toxicity) 
                          & !is.na(aquatic_chronic))

model_2<-polr(factor~ year_new + 
                toxicity+ 
                toxicity_worse+
                cancer_epa +
                aquatic_chronic +
                aquatic_acute, data=data_safe_cfrfull)

model_full <-polr(factor~ year_new + 
                    toxicity+ 
                    toxicity_worse+
                    cancer_epa +
                    aquatic_chronic +
                    aquatic_acute+
                    fruit_veg +
                    one_of_3_primary_acreage_crops_US + 
                    one_of_3_primary_acreage_crops_US*year_new
                  , data= data_safe_cfrfull)


###models with robust standard errors###
vcov_1 <- vcovCL(model_1, as.character(data_safe_cfr$pesticide))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_2, as.character(data_safe_cfrfull$pesticide))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_full, as.character(data_safe_cfrfull$pesticide))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_1,model_2,model_full, se=list(new.ses1, new.ses2, new.sesfull))


###ROBUSTNESS CHECK 2 (Table 4) - NO REVOKED####
datacfr$factorn<-datacfr$factor
datacfr$factorn[datacfr$factorn==4 ] <- NA
datacfr$factorn<-as.factor(datacfr$factorn)

######models for paper#####
datacfr_norev<-subset(datacfr,!is.na(factorn) & !is.na(year_new))
model_1 <- polr(factorn~ year_new, data=datacfr_norev)

datacfr_norev_full<-subset(datacfr_norev,!is.na(toxicity_worse) 
                           & !is.na(cancer_epa) 
                           & !is.na(toxicity) 
                           & !is.na(aquatic_chronic))

model_2<-polr(factorn~ year_new + 
                toxicity+ 
                toxicity_worse+
                cancer_epa +
                aquatic_chronic +
                aquatic_acute, data=datacfr_norev_full)

model_full <-polr(factorn~ year_new + 
                    toxicity+ 
                    toxicity_worse+
                    cancer_epa +
                    aquatic_chronic +
                    aquatic_acute+
                    fruit_veg +
                    one_of_3_primary_acreage_crops_US + 
                    one_of_3_primary_acreage_crops_US*year_new
                  , data= datacfr_norev_full)


###models with robust standard errors###
vcov_1 <- vcovCL(model_1, as.character(datacfr_norev$pesticide))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_2, as.character(datacfr_norev_full$pesticide))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_full, as.character(datacfr_norev_full$pesticide))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_1,model_2,model_full, se=list(new.ses1, new.ses2, new.sesfull))


####ROBUSTNESS CHECK WITH EU CANCER (APPENDIX TABLE 1)###
model_eu1 <- polr(factor~ year_new, data=datacfr1)

model_eu2<-polr(factor~ year_new + 
                  toxicity+ 
                  toxicity_worse+
                  cancer_echa +
                  aquatic_chronic +
                  aquatic_acute, data=datacfrfull)

model_eufull <-polr(factor~ year_new + 
                      toxicity+ 
                      toxicity_worse+
                      cancer_echa+ 
                      aquatic_chronic +
                      aquatic_acute+
                      fruit_veg +
                      one_of_3_primary_acreage_crops_US + 
                      one_of_3_primary_acreage_crops_US*year_new
                    , data= datacfrfull)


###models with robust standard errors###
vcov_1 <- vcovCL(model_eu1, as.character(datacfr1$pesticide))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_eu2, as.character(datacfrfull$pesticide))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_eufull, as.character(datacfrfull$pesticide))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_eu1,model_eu2,model_eufull, se=list(new.ses1, new.ses2, new.sesfull))

####ROBUSTNESS CHECK WITH HIGHER ORDER CLUSTERING (APPENDIX TABLE 2)###
model_1 <- polr(factor~ year_new, data=datacfr1,Hess=TRUE)

model_2<-polr(factor~ year_new + 
                toxicity+ 
                toxicity_worse+
                cancer_epa +
                aquatic_chronic +
                aquatic_acute, data=datacfrfull, Hess=TRUE)

model_full <-polr(factor~ year_new +
                    toxicity+ 
                    toxicity_worse+
                    cancer_epa +
                    aquatic_chronic +
                    aquatic_acute+
                    fruit_veg +
                    one_of_3_primary_acreage_crops_US + 
                    one_of_3_primary_acreage_crops_US*year_new
                  , Hess=TRUE,data= datacfrfull)


###models with robust standard errors###
vcov_1 <- vcovCL(model_1, as.character(datacfr1$tox_grouping))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_2, as.character(datacfrfull$tox_grouping))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_full, as.character(datacfrfull$tox_grouping))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_1,model_2,model_full, se=list(new.ses1, new.ses2, new.sesfull))


###ROBUSTNESS CHECK WITH AGE DUMMY (APPENDIX TABLE 3)###
datacfr$agedummies <-datacfr$year_new
datacfr$agedummies[datacfr$agedummies<30] <- 0
datacfr$agedummies[datacfr$agedummies>=30] <- 1
datacfr$agedummies<-as.factor(datacfr$agedummies)

##models for appendix

model_1_dum <- polr(factor~ agedummies, data=datacfr)

datacfrfull_dum<-subset(datacfr,!is.na(toxicity_worse) 
                        & !is.na(cancer_epa) 
                        & !is.na(toxicity) 
                        & !is.na(aquatic_chronic))

model_2_dum<-polr(factor~ agedummies + 
                    toxicity+ 
                    toxicity_worse+
                    cancer_epa +
                    aquatic_chronic +
                    aquatic_acute, data=datacfrfull_dum)

model_full_dum <-polr(factor~ agedummies + 
                        toxicity+ 
                        toxicity_worse+
                        cancer_epa +
                        aquatic_chronic +
                        aquatic_acute+
                        fruit_veg +
                        one_of_3_primary_acreage_crops_US 
                      , data= datacfrfull_dum)


###models with robust standard errors###
vcov_1 <- vcovCL(model_1_dum, as.character(datacfr$pesticide))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_2_dum, as.character(datacfrfull_dum$pesticide))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_full_dum, as.character(datacfrfull_dum$pesticide))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_1_dum,model_2_dum,model_full_dum, se=list(new.ses1, new.ses2, new.sesfull))


###ROBUSTNESS CHECK WITH 3 AGE BUCKETS (APPENDIX TABLE 4)###
datacfr$agedummies <-datacfr$year_new
datacfr$agedummies[datacfr$agedummies<30] <- 0
datacfr$agedummies[datacfr$agedummies>=30&datacfr$agedummies<46] <- 1
datacfr$agedummies[datacfr$agedummies>=46] <- 2
datacfr$agedummies<-as.factor(datacfr$agedummies)

##models for appendix

model_1_dum <- polr(factor~ agedummies, data=datacfr)

datacfrfull_dum<-subset(datacfr,!is.na(toxicity_worse) 
                        & !is.na(cancer_epa) 
                        & !is.na(toxicity) 
                        & !is.na(aquatic_chronic))

model_2_dum<-polr(factor~ agedummies + 
                    toxicity+ 
                    toxicity_worse+
                    cancer_epa +
                    aquatic_chronic +
                    aquatic_acute, data=datacfrfull_dum)

model_full_dum <-polr(factor~ agedummies + 
                        toxicity+ 
                        toxicity_worse+
                        cancer_epa +
                        aquatic_chronic +
                        aquatic_acute+
                        fruit_veg +
                        one_of_3_primary_acreage_crops_US 
                      , data= datacfrfull_dum)


###models with robust standard errors###
vcov_1 <- vcovCL(model_1_dum, as.character(datacfr$pesticide))
new.ses1 <- sqrt(diag(vcov_1))
vcov_2 <- vcovCL(model_2_dum, as.character(datacfrfull_dum$pesticide))
new.ses2 <- sqrt(diag(vcov_2))
vcov.full <- vcovCL(model_full_dum, as.character(datacfrfull_dum$pesticide))
new.sesfull <- sqrt(diag(vcov.full))

stargazer(model_1_dum,model_2_dum,model_full_dum, se=list(new.ses1, new.ses2, new.sesfull))


####PETITION DATA REPLICATION###
data<-read.csv("petition_replication_data.csv")

##SUMMARY STATS (TABLE 5)
#Company petitions
co_petitions <- subset(data,petitioner!="IR4"& petitioner!="ciel" & petitioner!="y-tex" & petitioner!="pace" & petitioner!= "kim-ci and siemer")
co_petitions$change[co_petitions$change>0] <- 0 #change is requested change
co_petitions$change[co_petitions$change<0] <- 1
co_petitions$change<-as.factor(co_petitions$change)
summary(co_petitions$change)/length(co_petitions$change)

#IR4 petitions
ir4_petitions <- subset(data,petitioner=="IR4")
ir4_petitions$change[ir4_petitions$change>0] <- 0
ir4_petitions$change[ir4_petitions$change<0] <- 1
ir4_petitions$change<-as.factor(ir4_petitions$change)
summary(ir4_petitions$change)/length(ir4_petitions$change)

#pesticide age at time of petition for IR4 and company petitions
co_petitions$age_co <- co_petitions$year_published_fr - 
  co_petitions$year_registered
ir4_petitions$age_ir4 <- ir4_petitions$year_published_fr - 
  ir4_petitions$year_registered

##MODEL 1 COMPANIES TABLE 6###
co_model<-glm(co_petitions$change ~ co_petitions$age_co, na.action=na.exclude,family="binomial")
###Robust errors##
vcov.matcfr1 <- vcovCL(co_model, as.character(co_petitions$pesticide))
new.ses1 <- sqrt(diag(vcov.matcfr1))
stargazer(co_model,se= new.ses1)

##MODEL 2 COMPANIES TABLE 6###
co_petitions$di_age<-co_petitions$age_co
co_petitions$di_age[co_petitions$di_age<=20]<-0
co_petitions$di_age[co_petitions$di_age>20]<-1
di_co_model<-glm(co_petitions$change ~ co_petitions$di_age, na.action=na.exclude,family="binomial")
###Robust errors##
vcov.matcfr1 <- vcovCL(di_co_model, as.character(co_petitions$pesticide))
new.ses1 <- sqrt(diag(vcov.matcfr1))
stargazer(di_co_model,se= new.ses1)


##MODEL 1 IR4 TABLE 6###
ir4_petitions <-subset(ir4_petitions,!is.na(age_ir4) & !is.na(change))
model<-glm(change ~ age_ir4, data= ir4_petitions, na.action=na.exclude, family="binomial")
vcov.matcfr1 <- vcovCL(model, as.character(ir4_petitions $pesticide))
new.ses1 <- sqrt(diag(vcov.matcfr1))
stargazer(model, se=new.ses1)

##MODEL 2 IR4 TABLE 6###
ir4_petitions$di_age<-ir4_petitions$age_ir4
ir4_petitions$di_age[ir4_petitions$di_age<=20]<-0
ir4_petitions$di_age[ir4_petitions$di_age>20]<-1
di_ir4_model<-glm(ir4_petitions$change ~ ir4_petitions$di_age, na.action=na.exclude,family="binomial")
vcov.matcfr1 <- vcovCL(di_ir4_model, as.character(ir4_petitions$pesticide))
new.ses1 <- sqrt(diag(vcov.matcfr1))
stargazer(di_ir4_model, se=new.ses1)


##COMPANY PETITIONS WITH DICHOTOMOUS AGE AND VARIOUS CUT POINTS (APPENDIX TABLE 5)###
#model 2
co_petitions$di_age<-co_petitions$age_co
co_petitions$di_age[co_petitions$di_age<=19]<-0
co_petitions$di_age[co_petitions$di_age>19]<-1
di_co_model<-glm(co_petitions$change ~ co_petitions$di_age, na.action=na.exclude,family="binomial")
###Robust errors##
vcov.matcfr1 <- vcovCL(di_co_model, as.character(co_petitions$pesticide))
new.ses1 <- sqrt(diag(vcov.matcfr1))
stargazer(di_co_model,se= new.ses1)


#model 3
co_petitions$di_age<-co_petitions$age_co
co_petitions$di_age[co_petitions$di_age<=18]<-0
co_petitions$di_age[co_petitions$di_age>18]<-1
di_co_model<-glm(co_petitions$change ~ co_petitions$di_age, na.action=na.exclude,family="binomial")
###Robust errors##
vcov.matcfr1 <- vcovCL(di_co_model, as.character(co_petitions$pesticide))
new.ses1 <- sqrt(diag(vcov.matcfr1))
stargazer(di_co_model,se= new.ses1)
